///@file Vector.cpp
///@brief Implementation of the class representing a Vector
///@author Arnaud Duval
///@version 0.0
///@date 2010/11/10



#include <cmath>
#include <cstring>
#include "Vector.h"


///@brief Default constructor
Vector::Vector()
{
	size = 0;
}

///@brief Constructor precising the size of the vector
///@param n Size of the vector
///@param val Default value of vector components
Vector::Vector(unsigned int n, double val)
{
	assert(n > 0);
	size = n;
	a = new double[size];
	///@brief Fill vector with default value
	for(unsigned int i = 0 ; i < size ; i++)
		a[i] = val;
}

///@brief Copy constructor
///@param v vector to copy
Vector::Vector(const Vector &v)
{
	size = v.size;
	a = new double[size];
	for(unsigned int i = 0 ; i < size ; i++)
		a[i] = v.a[i];
}

///@brief Copy constructor for temporary object
///@param v Pointer to vector to copy
Vector::Vector(Vector* v)
{
	size = (*v).size;
	a = new double[size];
	for(unsigned int i = 0 ; i < size ; i++)
		a[i] = (*v).a[i];
}

///@brief Destructor
Vector::~Vector()
{
	if(size != 0)
		delete[] a;
}

///@brief Overload of = operator
///@param v right hand member vector
Vector& Vector::operator=(const Vector &v)
{
	if(size)
		delete[] a;
	size = v.size;
	a = new double[size];
	for(unsigned int i = 0 ; i < size ; i++)
		a[i] = v.a[i];
	return (*this);
}


///@brief Overload of = operator for temporary object
///@param v Right hand operand
Vector& Vector::operator=(Vector* v)
{
	size = (*v).size;
	a = new double[size];
	for(unsigned int i = 0 ; i < size ; i++)
		a[i] = (*v).a[i];
	return (*this);
}

///@brief Overload of &= operator : copy without size check
///@param v Vector object to copy
Vector& Vector::operator&=(const Vector & v)
{
	if(!v.size)
	{
		if(size)
			delete[] a;
		size = 0;
	}
	else
	{
		if(v.size != size)
		{
			if(size)
				delete[] a;
			a = new double[size = v.size];
		}
		memcpy(a, v.a, v.size*sizeof(double));
	}
	return (*this);
}

///@brief Overload of [] operator : direct access to vector components
///@param i index of returned component
double& Vector::operator[](unsigned int i) const
{
	assert((i >=1)&&(i <= size));
	return a[i-1];
}

///@brief Overload of += operator
///@param v right hand operand
void Vector::operator += (const Vector& v)
{
	assert (size==v.size);
	for (unsigned int i=0; i<size; ++i)
		a[i] += v.a[i];
}
///@brief Overload of -= operator
///@param v right hand operand
void Vector::operator -= (const Vector& v)
{
	assert (size==v.size);
	for (unsigned int i=0; i<size; ++i)
		a[i] -= v.a[i];
}

///@brief Overload of /= operator
///@param d Double precision value. Only exact 0 value is prevented
void Vector::operator /= (const double& d)
{
	assert(d != 0.);
	for(unsigned int i=0 ; i <size ; ++i)
		a[i] /= d;
}


///@brief Compute euclidian norm of a vector
double Vector::Euclide_Norm()
{
	double sum = 0.;
	for(unsigned int i = 0 ; i < size ; i++)
		sum += a[i]*a[i];	
	sum = sqrt(sum);
	return sum;
}

///@brief Return dimension of vector	
unsigned int Vector::Dim() const
{
	return size;	
}

///@brief Resizes the vector
///@param n New size of the vector
void Vector::resize(unsigned int n)
{
	assert(n > 0);
	if(size > 0)
		delete[] a;
	size = n;
	a = new double[n];
}
	
///@brief Swap 2 components
///@param i1 Index of first component
///@param i2 Index of scond component
void Vector::swap(unsigned int i1, unsigned int i2)
{
	assert ((i1 > 0)&&(i2>0));
	assert ((i1 <= size)&&(i2 <= size));
	double temp = a[i1 - 1];
	a[i1 - 1] = a [i2 - 1];
	a[i2 - 1] = temp;
}

///@brief Add a subvector into current vector
///@param v Subvector to add
///@param id Index of component where the first component of subvvector will be placed
void Vector::AddSubVector(Vector v, unsigned int id)
{
	assert(v.Dim()+id-1 <= size);
	for(unsigned int i = 0 ; i < v.Dim() ; i++)
	{
		a[id+i-1] += v[i+1];
	}
}

///@brief Fill vector with a Sobol sequence on [0,1]
void Vector::FillSobol()
{
	assert(size != 0);
	for(unsigned int i = 1 ; i <= size ; i++)
	{
		a[i-1] = 0.;
		
		int n=i;
		int count = 1;
		while(n!=0)
		{
			a[i-1] += ((double)(n%2))/(pow(2., count));
			n = n/2;
			count ++;
		}
	}
}

///@brief Sort the vector values with bubble sort alogorithm
void Vector::BubbleSort()
{
	assert(size > 0);
	bool permut = true;
	while(permut)
	{
		permut = false;
		for(unsigned int i = 1 ; i < size ; i++)
		{
			if(a[i-1] > a[i])
			{
				double temp = a[i];
				a[i] = a[i-1];
				a[i-1] = temp;
				permut = true;
			}
		}
	}
}

///@brief Generate permutations
///@param b map of indexes
///@param reset boolean to set the indexes map for the first permutation
Vector Vector::GetNextPermutation(int *b, bool reset)
{
	assert(size != 0);

	if(reset)
	{
		for(unsigned int j = 0 ; j < size ; j++)
			b[j] = j;
	}

	//Test if maximum number of permutation has been reached
	bool reached = true;
	for(unsigned int j = 0 ; j < size ; j++)
	{
		if(b[j] != size - j - 1)
			reached = false;
	}
	if(reached)
		return *(this);

	int temp;
	//Find largest index j with b[j] < b[j+1]
	int j = size - 2;
	while(b[j] > b[j+1])
		j--;

	//Find next index k such that b[k] is smallest integer greater than b[j] to the right of a[j]
	int k = size - 1;
	while(b[j] > b[k])
		k--;

	//Interchange b[j] and b[k]
	temp = b[k];
	b[k] = b[j];
	b[j] = temp;

	//Put tail end of permutation after jth position in increasing order
	int r = size - 1;
	int s = j + 1;

	while(r > s)
	{
		temp = b[s];
		b[s] = b[r];
		b[r] = temp;
		r--;
		s++;
	}

	Vector res(size);
	for(unsigned int j = 0 ; j < size ; j++)
		res.a[j] = a[b[j]];

	return res;
}

///@brief Return maximum value of Vector terms
double Vector::Max()
{
	assert(size >= 1);
	double res = a[0];
	for(unsigned int i = 1 ; i < size ; i++)
	{
		if(a[i] > res)
			res = a[i];
	}
	return res;
}
///@brief Return the sum of Vector terms
double Vector::Sum()
{
	assert(size >= 1);
	double res = 0.;
	for(unsigned int i = 0 ; i < size ; i++)
		res += a[i];
	return res;
}

///@brief Return average value of Vector components
double Vector::Average()
{
	assert(size >= 1);
	return Sum()/size;
}

///@brief Set all components to zero
void Vector::Nullify()
{
    assert(size >= 1);
    for(unsigned int i = 0 ; i < size ; i++)
        a[i] = 0.;
}

///@brief Overload of + operator
///@param v1 left hand member vector
///@param v2 right hand member vector
Vector operator + (const Vector& v1, const Vector& v2)
{
    assert(v1.size == v2.size);
    Vector res(v1.size);
    for(unsigned int i = 0 ; i < v1.size ; i++)
        res.a[i] = v1.a[i] + v2.a[i];
    return res;
}

///@brief Overload of - operator
///@param v1 left hand member vector
///@param v2 right hand member vector
Vector operator - (const Vector& v1, const Vector& v2)
{
	assert(v1.size == v2.size);
	Vector res(v1.size);
	for(unsigned int i = 0 ; i < v1.size ; i++)
		res.a[i] = v1.a[i] - v2.a[i];
	return res;	
}


///@brief Overload of | operator : scalar product
///@brief Computr \f$ v1_{i} v2_{i} \f$
///@param v1 left hand member vector
///@param v2 right hand member vector
double operator | (const Vector& v1, const Vector& v2)
{
	assert(v1.size == v2.size);
	double res = 0.;
	for(unsigned int i = 0 ; i < v1.size ; i++)
		res += v1.a[i]*v2.a[i];
	return res;	
}

///@brief Overload of * operator : product vector-scalar
///@brief Compute \f$ a_{i} = v_{i} \cdot s \f$
///@param v left hand member vector
///@param s right hand member scalar
Vector operator * (const Vector& v, const double& s)
{
	Vector res(v.size);
	for(unsigned int i = 0 ; i < v.size ; i++)
		res.a[i] = v.a[i]*s;
	return res;
}

///@brief Overload of * operator : product scalar-vector
///@brief Compute \f$ a_{i} = s \cdot v_{i} \f$
///@param s left hand member scalar
///@param v right hand member vector
Vector operator * (const double& s, const Vector& v)
{
	Vector res(v.size);
	for(unsigned int i = 0 ; i < v.size ; i++)
		res.a[i] = s*v.a[i];
	return res;	
}

///@brief Overload of / operator : division of a vector by a scalar
///@brief Compute \f$ a_{i} = frac{v_{i}}{s}
///@param v left operand : vector
///@param s right operand : scalar
///@warning This operation is not protected against divisions by very low value number. Use with caution
Vector operator / (const Vector& v, const double& s)
{
	Vector res(v.size);
	for(unsigned int i = 0 ; i < v.size ; i++)
		res.a[i] = v.a[i] / s;
	return res;
}

///@brief Overload of && operator : happend a vector to another one
///@param v1 first vector
///@param v2 second vector
Vector operator && (const Vector& v1, const Vector& v2)
{
	Vector res(v1.size + v2.size);
	for(unsigned int i = 0 ; i < v1.size ; i++)
	{
		res.a[i] = v1.a[i];
	}
	for(unsigned int i = 0 ; i < v2.size ; i++)
	{
		res.a[i+v1.size] = v2.a[i];	
	}
	
	return res;
}

///@brief Overload of << operator
///@param f out stream
///@param v Vector object to put in the out stream
std::ostream& operator<<(std::ostream &f, const Vector& v)
{
	for (unsigned int i = 0 ; i < v.size ; i++)
	{
		f << v.a[i] << "\t";
	}
	return f;
}

///@brief Overload of >> operator
///@param f in stream
///@param v Vector object to put in the in stream
std::istream& operator>>(std::istream &f, const Vector&v)
{
	for(unsigned int i = 0 ; i < v.size ; i++)
	{
		f >> v.a[i] ;
	}
	return f;
}

